Continuum limit of amorphous elastic bodies: 
A finite-size study of low frequency harmonic vibrations 
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The approach of the elastic continuum limit in small amorphous bodies formed by weakly poly- 
disperse Lennard- Jones beads is investigated in a systematic finite-size study. We show that classical 
continuum elasticity breaks down when the wavelength of the sollicitation is smaller than a charac- 
teristic length of approximately 30 molecular sizes. Due to this surprisingly large effect ensembles 
containing up to N = 40, 000 particles have been required in two dimensions to yield a convincing 
match with the classical continuum predictions for the eigenfrequency spectrum of disk-shaped ag- 
gregates and periodic bulk systems. The existence of an effective length scale £ is confirmed by the 
analysis of the (non-gaussian) noisy part of the low frequency vibrational eigenmodes. Moreover, 
we relate it to the non-affine part of the displacement fields under imposed elongation and shear. 
Similar correlations (vortices) are indeed observed on distances up to £ m 30 particle sizes. 

Pacs: 72.80.Ng, 65.60.+a, 61.46. +w 
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I. INTRODUCTION. 

Determining the vibration frequencies and the asso- 
ciated displacement fields of solid bodies with variaus. 
shapes is a well studied area of continuum mechanicsErEl 
with applications in fields as different as planetary science 
and nuclear physics. The increasing development of ma- 
terials containing nanometric size structures leads one to 
question the limits of applicability of classical continuum 
elasticity theory, which is in principle valid only on leni 
scales much larger than the interatomic distances! 
This question is relevant from an experimental viewpoint, 
since mechanical properties are inferred from spectro- 
scopic measurements systematically interpreted within 
the framework of continuum elasticityETu. As .increas- 
ingly smaller length scales are now investigated^, direct 
verification of this assumption is highly warranted. 

For macroscopic systems, on the other side, it is well 
known that the vibrational density of states in amor- 
phous glassy materials deviates from the classical spec- 
trum at the so called ".-Bason peak" frequency which is 
in the Terahertz rangeuTL3. i-3jhe nature of the "Boson 
peak" is highly controversialllJlla. However, this experi- 
mental fact suggests that continuum theory is inappro- 
priate at small length scales where the disorder of the 
amorphous system may become relevantLVj. Obviously, 
one may ask if there is a finite length scale below which 
the classical mechanical approach becomes inappropri- 
ate, and-Hihat the microscopic features are which deter- 
mine itc tfu. Elaborating further the brief presentation 
given ir£3, we show by means of a simple generic simula- 
tion model that, indeed, a relatively large characteristic 
length exists, and, second, that it envolves collective par- 
ticle rearrangements. 



The above questions are more generally related to the 
propagatioii-|Of waves in disordered materialstll, and con- 
cern^oamnJ and emulsionscil as well as granular mate- 
rialsta~E3, when they are submitted to small amplitude 
vibrations. As for these systems of current interest, the 
existence of an elastic limit is still a matter of debateE3, 
we believe that the detailed characterization of strongly 
heterogeneous elastic benchmark systems with a well de- 
fined continuum limit approach is crucial. 

In this paper, we investigate the existence of a contin- 
uum limit in the vibrational modes of two-dimensional 
amorphous nanometric Lennard- Jones materials. The 
objects we consider are either disk-shaped clusters of di- 
ameter 2R, as the one shown on the left side of Fig. [l], or 
bulk like systems without surfaces contained in a square 
of side L with periodic boundary conditions (Fig. [j](b)). 
Technically, the systems are formed by carefully quench- 
ing a slightly polydisperse liquid of spherical particles in- 
teracting via simple Lennard- Jones (LJ) pair potentials 
into the nearest energy minimum. Due to the polydisper- 
sity the resulting structures are isotropic and amorphous, 
i.e. exhibit no long range crystalline order. The force net- 
work (Fig. [I]) appears to be strongly varying with weak 
and tensile zones (red) embedded within a rigid repulsive 
skeleton (black). These "force chains" are very similar to 
those found in_cohesionless granular media without at- 
tractive forcesEj. This feature may be added to the list 
of similarities which have been noticed pbetween granu- 
lar and amorphous (glassy) materialsE3e.3. As the force 
network is strongly inhomogeneous, the relevance of the 
quenched stresses is a natural questionQ. 

We investigate the vibrational modes of these objects 
using atomic level simulations. All particle coordinates 
and interparticle forces are exactly known here, and it is 
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possible to calculate the vibration frequencies around an 
equilibrium position, by— exact diagonalization of the so 
called dynamical matrfxnJ expressible in terms of the first 
and second derivatives of the interparticle interaction po- 
tentials. We have carried out a systematic comparison of 
these eigenfrequencies w(p) (p being an index increas- 
ing with frequency) obtained numerically with those pre- 
dicted by continuum elasticity-for two dimensional ob- 
jects of increasingly large sizesEJ. We concentrate on the 
lowest end of the vibrational spectrum, since this is the 
part that corresponds to the largest wavelengths for the 
vibrations, and should reach first the continuum limit. 
These frequencies are also those which arej-probed in low 
frequency Raman scattering experimentsO, in order to 
determine the typical size of nanoparticles. 

The key result of this paper is to show the existence 
of a characteristic wavelength (thus a characteristic size) 
beyond which the classical continuum limit is valid, but 
below which it is erroneous. Moreover, we show the ex- 
istence of rotational structures (vortices) of similar sizes, 
when the system is submitted to simple mechanical sol- 
licitations (traction and shear). The size of these vortices 
is relatively large (« 30 average interatomic distances). 
We discuss the relation between the sizes of the vortices 
and the limit of applicability of the classical continuum 
theory by computing the elastic moduli, by studying the 
symmetry of the nanoscale stress tensor, and by identi- 
fying the low frequency vibrational eigenmodes. 

Our paper is arranged as follows: In Sec. [n] we sum- 
marize some basic relations and results of classical con- 
tinuum theory for two dimensional clastic bodies. Sim- 
ulation techniques, sample parameters and preparation 
protocols of our model amorphous systems are explained 
in Sec. HI and simple system properties are discussed. In 



Sec. IV we analyse histograms and spatial correlations of 
the quenched forces and of the stiffness of the bonds. A 
weak enhancement of the rigid skeleton is demonstrated 
for small systems. The next two sections contain the key 
results of this paper. In Sec. [v] we discuss the mechani- 
cal properties of a periodic bulk system under elongation 
and shear, we compute the elastic moduli and we char- 
acterize the non-affine displacement field generated. The 
eigenvalues and eigenvectors and their departures from 
the continuum prediction are analyzed in Sec. VI. We 
conclude with a summary of our results in Sec. VII. 



II. THEORY 



Before presenting our computational results we give 
here a short synopsis of assumptions and results of clas- 
sical continuum elasticity in two dimensions, relevant for 
the subsequent discussion of our numerical results. 



A. Continuum description of an isotropic elastic 
body in 2D. 

The main assumption of classical elasticity^, is that 
the system can be entirely described by a unique vectorial 
field: the displacement field u(r) , describing the displace- 
ment of a volume element from its equilibrium position r. 
As we are interested in systems at zero temperature, en- 
ergy and fxee energy are obviously identical. The Landau 
expansionB of the energy 8E per unit volume is expressed 
in terms of u and its first derivatives. Due to translational 
and rotational invarianceEE of 5E, it depends only up to 
second order on the symmetric part of grad u, the (lin- 
earized) strain tensor e Q /3 = 1/2 {dua/dxp + dup/dx a ) 
with x\ = x and xi = y for the coordinates. Moreover in 
linear elasticity (that is in the lowest order in the Landau 
expansion), for isotropic and homogeneous systems only 
two Landau parameters are required, thus: 
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5E=^ (e xx + e yy f + n (e 2 xx + e 2 yy + 2e 2 xy ) (1) 



where we have defined the two phenomenological Lam 
coefficients A and [i. fj, is the so called shear modulus, 
and A takes part in the compressibility k = A + fi. Be- 
low, we will discuss two methods to determine them in a 
computer experiment. The stress tensor can be defined 
as the conjugate variable of the strain tensor, 



cr Q/3 = d5E/de a(L 



(2) 



In this definition a a p it obviously a symmetric tensor. 
We will see later that a microscopically constructed stress 
tensor can, however, violate this symmetry condition at 
small length scales (section |ll C| ). 

For the remaining three independent stress tensor el- 
ements, one gets from Eqns (|l|) and ( |^) the standard 
Hooke 's relations: 
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(A + 2fx)e xx + \e y 



xx — v v ' ^h^J^xx i / ^yy 
(T yy = (A + 2/j,)e yy + \e f , 
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We use these formulae in Sec. [y] to measure directly the 
Lame coefficients from the forces generated in a periodic 
box under external strain. It is easy to work out from 
Eq. ( |^) that in 2 dimensions the Poisson ratio is given 
by v = XJJX + 2fi) with an upper bound 1 (and not 1/2 
as in 3D)0. 

As it is well knownEm, the equations of motion 
pda a p/dx a — u a together with the constitute Eqns. ( |^) 
correspond to wave equations with boundary conditions 
depending on the problem of interest: for the periodic 
box problem the solutions must have the same period- 
icity, and the freely floating disk-shaped aggregate re- 
quires vanishing lateral and radial stresses on its surface. 
The velocities of transverse and longitudinal waves which 
solve these wave equations are given in terms of the Lame 
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coefficients, i.e. one has c\ =nAV 'p an( l c l = + ^p)/p 
where p is the particle densitytl. We remind that cl > ct 
and that transverse modes thus correspond generally to 
smaller cigcnfrequencies. 

The solutions for the periodic box case are of course the 
plane waves with a wavevector quantified by the bouad- 
ary conditions, (k x ,k y ) — ^-(n,m), with wavelengthc^l 

and with dimcnsionless frequency 

with two quantum numbers n, m = 0, f,.... The run- 
ning index p increases with frequency. In the continuous 
case, the dispersion relation is linear, and the frequency 
is straightforward. Hence, eigenfrequencies are charac- 
terized by a pair of different integers. They are eight- 
fold degenerated if n ^ m ^ and four-fold in all other 
cases. The associated plane waves travel in two opposite 
and orthogonal directions. 

The situation for disk-shaped objects is somewhat 
more complexa, with again two quantum numbers n and 
k characterizing the cigcnmodes. The quantum number 
k is associated with the angular dependency of the dis- 
placement field (and is due to the 2tt periodicity), and the 
number n to its radial dependency. The eigenfrequencies 
are obtained by solving the non linear dispersion relation. 
They are of the form 

^(p)-(^f) 2 = U{u) (6) 

where v is the Poisson ratio given above. The eigenvec- 
tors are related to Bessel functions^ which may be ap- 
proximated locally by plane waves if X(p) <C R. Vibra- 
tional modes of disks are either non degenerate (for axi- 
ally symmetric modes) or have two-fold degeneracy (for 
all other modes). Indeed, for k > 0, for every solution 
whose amplitude is oc e lke , one finds a second solution, 
orthogonal to the first one, with amplitude oc e lk ( s +^) 
where A8 = ir/2k. Every additional solution with the 
same A: is a linear combination of these two vectors. For 
axially symmetric modes (fc = 0), the above argument 
does not apply; any additional solution found by turning 
the coordinate system is identical to the first one. 

We finally stress the obvious: degenerate eigenvalues 
are inherent to the continuum treatment of highly sym- 
metric systems. The failure to observe them indicates 
either the lifting of the relevant symmetry or the break- 
down of continuum theory. 



B. Pair potential systems with central forces. 

In a simulation one has the advantage to know all the 
individual contributions to the total energy. The situa- 
tion is particularly simple if one has to deal with inter- 
particle pair potentials U (ry ) (ry being the interparticle 
distance) such as the LJ potential we use, and if we stay 
at zero temperature. In this case, the difference of the 
total potential energy 

N-l N 
8=1 j>i 

due to a displacement field u can be written to second 
order as a Hessian form 8E p = in* ■ M ■ u in terms of 
the (2N) x (27V) dynamical matrix M whose elements are 
given by 

- r^jMiajp = i (8 a p - n a n/3) + r?jCijn a np (8) 

n being the unit vector of the bond (for simplicity, we 
do not indicate the dependence of n on the particle 
indices i and j), tij = dU(rij)/drij the tension and 
ay = d 2 U(rij)/drfj the stiffness of the bond between 
two interacting beads i and j. The first one is related to 
the stresses "quenched" in the bulk, but as we will see 
later it is in general small in comparison with the sec- 
ond one. As in the present study the mass m of each 
monomer is set equal to unity (even though the particle 
diameters are polydisperse) the Euler-Lagrange equation 
M ■ u + mu = is solved directly by eigenfrequencies and 
displacement fields obtained simply by diagonalisation of 
the dynamical matrix. 

Assuming constant deformations (purely affine dis- 
placement field) under elongation and shear, one may 
use the dynamical-^iiatrix to calculate the Lame coeffi- 
cients. Comparing^ the Landau expression Eq. ( |l|) and 
the microscopic energy information in Eq. ( we get 

Xa = I E V^ n >l + r^Cij ( « + n A y )/2 - 2n 2 x n 2 y )] 

ij 

l l a = ~| ^ t r ij c ij n x n y (9) 
ij 

with A being the total surface. The sums run over all 
pairs of particles. As an affine displacement field is as- 
sumed to be valid down to atomic distances the coeffi- 
cients have been assigned an index a to distinguish them 
from the true macroscopic Lame coefficients that will be 
calculated later. It is instructive to consider the simple 
one dimensional analogy, i.e. a long string of connected 
springs with force constants q, to appreciate the approx- 
imation made in the above equations. We recall that 
the effective macroscopic Young modulus is A = l/(l/c;) 
rather than A a = (cj). This is due to the fact that the 
tension along the chain is constant but not the displace- 
ment field. Hence, Eqns. ( ^|) arc in general only good 
approximations if the dispersion of the effective spring 
constants is small and an external macroscopic strain 
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is locally transmitted in an affine way. (A is then in- 
deed given by the mean coupling constant with a neg- 
ative quadratic dispersion correction in leading order.) 
We will test this assumption in Sec. and estimate the 
length below which the affinity approximation becomes 
problematic. 



C. Definition and Measurement of stress tensor. 



In the elastic continuum theory (section [I A), the 
stress tensor is usually introduced by considering the 
forces per unit length exerted between adjacent volume 
elements. The total energy 6E thus corresponds to the 
work of internal forces. At a microscopic level, a natural 
way of extending this macroscopic approach is to com- 
pute a a /3 as the sum of forces per unit length (—tijn a /b) 
exerted in the a direction through the side perpendicular 
to the (3 direction of a small volume element of size b. A 
stress tensor defined in this way can be asymmetric. The 
usual, macroscopic proofirt] of the symmetry of a is based 
on the fact that intermolecular forces are short-ranged. 
Hence, one expects symmetry of the microscopic stress 
tensor (obtained from the 'macroscopic' definition) only 
for volume elements of size b much larger than the range 
of intermolecular forces. In section [V|, we will investigate 
how this macroscopic limit is reached. 

Note that our definition of the,stress tensor, does not 
correspond strictly to the usualoEa microscopic Kirkwood 
definition 



Vap{i) = —'Ejtijn a .np 



(10) 



which is necessarily symmetric. Although both quanti- 
ties yield the same macroscopic stress tensor, the defini- 
tion we are using is more appropriate to illustrate devia- 
tions from macroscopic behavior at small scales. Another 
method for illustrating these deviations can be found in 
Ref.tJ. In this work, the local breakdown of rotational 
invariance in a disordered system was shown to imply 
contribution of the asymmetric part of u to the local 
elastic energy. Defining o a p as 88E / d{du a / dj3) in place 
of Eq. ( ||), it can be shown that the contribution of the 
asymmetric part of a to the local energy is related to the 
work of the fluctuating part oLthe forces on the fluctu- 
ating part of u. In reference Ej, it was shown that this 
contribution vanishes only for very large systems. 



III. SAMPLE PREPARATION AND 
CHARACTERIZATION. 

Systems with two different boundary conditions (disks 
and periodic bulk systems) and three quench protocols 
have been simulated. In this section we discuss some 
technical points concerning the simulation methods and 
the sample preparation protocols and parameters. The 



details of the protocols and some properties of the final 
configurations are summarized in the Tables ffl, [n] and III . 

In the-, present study we use a shifted LJ 
potential for polydisperse particles with Ulj{t) — 
e ((r /r) 12 - 2(r /r) 6 ) + const for r < 2r and Ui,j{r) = 
otherwise, correctly regularized. We have written the 



potential in terms of the distance rg = (a^ 



0/2«/8 



where the LJ force vanishes. Natural LJ units are used, 
i.e. we set the energy parameter e = 1, the particle 
mass m = 1 and the mean diameter a — (a^) = 1. Note 
that while the particle mass is strictly monodisperse the 
particle diameters at are homogeneously distributed be- 
tween 0.8 and 1.2. This corresponds to a polydispersity 
index Sa/a ~ 0.12 which is sufficient to prevent large 
scale crystalline order. We did not attempt to make 
the particles even more polydisperse fearing demixing or 
systematic radial variation of particle sizes in the case of 
disk-shaped aggregates. 

The quench generally starts with molecular dynamics 
(MD) at some fixed .temperature using a simple velocity 
rescaling thermostattd. As integrator we use the veloc- 
ity Verlet algorithmEj with a time increment St = 0.001. 
The temperature remains constant over a fixed time in- 
terval of 1000 MD steps for each temperature step. This 
was sufficient to relax systems into a steady-state (obvi- 
ously, not necessarily the equilibrium) as monitored by 
pressure and system energy. Unfortunately, no more de- 
tailed characterization of the ageing as a function of the 
quench rate has been recorded. Following the initial MD 
sequence we quench the systems further down using a 
steepest descent (SD) method with a frictional force pro- 
portional to the total particle velocity and an additional 
damping term proportional to the relative velocity of in- 
teracting particles. We use a small value of order one 
for both friction constants and a simple leap-frog like ia- 
tegrator consistent with the velocity depended forcecJ. 
The particle motion is overdamped. The systems were 
cooled down with SD over a period of at least 1000 itinae. 
steps. Finally, the conjugate gradient method (CG)E3£-3 
was iterated till the configurations reach their local en- 
ergy minima. Double precision was used at every stage. 

An example for a disk-shaped aggregate has already 
been presented in fig. |l|(a). Two different protocols have 
been employed to generate such disks: 

• Protocol I: The starting point of the first protocol 
is a not too dense LJ liquid droplet of radius Ri 
and pi ~ 0.9 and temperature Ti = 0.1. The initial 
radius Ri is chosen such that the final mean radius 
R becomes not too different. We cool the systems 
with MD steps at T = 0.1,0.05,0.01,0.005 and 
0.001. As described above, we finally quench each 
system into its local minimum using a sequence of 
SD and CG steps. See Table | for details. 

• Protocol II: In the second case spherical disks of 
radius Ri were cut out of amorphous periodic sys- 
tems already prepared at T = (see Protocol III) 
and p — 0.925 and the disks are quenched as before 
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with SD and CG. No finite temperature MD step 
was included here. See Table [f| for details. The 
second protocol is much faster than the first one 
which might, however, mimic better the actually 
occurring aggregation process in real systems. 



The radius R — y2{Ti + h)/N of quenched cluster 
is obtained from the eigenvalues I\ > I2 of the inertia 
tensor of the cluster. The mean density p —< N/ttR 2 > 
obtained accordingly is a weakly decreasing function of 
the particle number as can be seen from the Tables | and 
O. This is expected and in qualitative agreement with 
the decreasing Laplace pressure P oc 1/R. Note that 
the pressure is very small for all the disks (Protocol I 
and II alike) and is not indicated. We have also tried 
to characterize the shape of the disks and have indicated 
the excentricity e = \J 1 — I2 / h ■ As can be seen, small 
aggregates are strongly elliptic, but the effect is much 
stronger for the first slow aggregation protocol. This is 
probably due to capillary waves formed at T = 0.1 which 
are subsequently frozen in. Obviously, ellipticity is one 
possible cause for lifting of the eigenfrequency degener- 
acy. This effect should, however, become small for larger 
clusters where e vanishes (Tab. || and ||). Moreover, we 
will see later that the mechanical properties of our aggre- 
gates do not depend significantly on the quench protocol. 

Periodic bulk systems, such as the one presented in 
Fig. ^(b), have been prepared following the third pro- 



tocol and are summarized in Table III. We started 
here with equilibrated liquids at T = 1 which were 
then subsequently cooled down with temperature steps 
at T = 0.5,0.1,0.05,0.01,0.005 and 0.001. Systems be- 
tween L = 7a and L — 208a and containing from N = 50 
up to N = 40, 000 particles have been sampled. 

For N = 10,000 particles we have systematically 
scanned over density varying the box size from L = 100 
to L = 110. This was done in order to find — for the 
given protocol — a working point density for which large 
bulk systems correspond to a near zero pressure state 
P(T = 0) w 0, thus the mechanical properties of the bulk 
systems correspond to the free-floating aggregates. Note 
that only one configuration has been sampled for this 
measurement sequence. As a sideline, we draw attention 
to the fact that systems with P < while mechanically 
stable in a periodic box are thermodynamically unstable 
at low temperature. If more time would be given to the 
systems to equilibrate as allowed by the protocol phase 
separation would occur. Indeed we find evidence for the 
formation of small holes for p < 0.85. In all other cases 
we find that density is perfectly homogeneous down to 
a scale typical of the interatomic distance and density 
fluctuations are easily excluded as a microscopic candi- 
date for explaining the slow continuum approach. Our 
systematic finite-size study, i.e. the variation of box sizes 
L, was performed at fixed density p = 0.925. In order 
to scale correctly the eigenfrequencies of the aggregates 
we have also recorded the Lame coefficients as a function 
of density. We have checked that the particle energies 



E, Lame factors A a and p a -, mean effective spring con- 
stants (rfjCij) (which are included in the different tables) 
do not depend significantly on the different symmetries 
and quench protocols, for a given density. It is reassuring 
that all the results we obtained are robust with regard to 
the variation of the quench protocol even though history 
dependence must obviously play a role for the details. 
Typically, we have generated 20 configurations for each L 
and in terms of CPU hours this was the most demanding 
part of this study. The simulations have been performed 
on a local workstation cluster over a period of one year. 

Interestingly, we notice in Table [nj that the pressure, 
the particle energy and the Lame coefficient A increase 
systematically in small boxes. This is shown specifically 
for A in Fig. [|(b) . These finite-size effects indicate corre- 
lations on a similar length even in larger systems. We will 
now turn to the characterization of these correlations, by 
first studying the dynamical matrix, and the quenched 
stresses. 



IV. CONTRIBUTIONS TO THE DYNAMICAL 
MATRIX: QUENCHED FORCES AND SPRING 
CONSTANTS 

In this section we discuss various distributions and cor- 
relation functions associated with the quenched particle 
positions, forces and spring constants which contribute to 
the dynamical matrix (Eq. (|8|)). We attempt a character- 
ization of the frozen in disorder visualized in Fig. [j] where 
the snapshots reveal strong spatially correlated fluctua- 
tions. Specifically, we ask if it is possible to extract a 
characteristic length scale solely from the distribution of 
weak and rigid regions which might be a candidate to 
explain the large crossover length scale £ mentioned in 
the Introduction. We start by giving some additional 
information concerning the snapshots, turn then to the 
histograms and discuss finally the fractal structure of the 
quenched forces. 

A timely justification for discussing quenched forces 
lies in the current interest of their role tt^ijajpular materi- 
als, foams and glassy colloidal systems&EirEl It has been 
suggested by S. Alexander!] and others that the quenched 
forces might contribute to the unusual mechanical and 
rheological properties in these systems. However, to put 
this immediately into perspective, quenched forces are 
unlikely candidates to rationalize alone the slow contin- 
uum approach discussed in Sec. 



VI 



basically, since their 
average contribution to the dynamical matrix is weak: 
(rfjCjj) 3> (fijtij). This inequality is revealed in Ta- 
ble HI. We come back to this in Sec. VI where we dis- 



cuss the contribution of the quenched forces to the eigen- 
modes. 

The second and more important point we want to make 
here is that for a given potential the first and the second 
derivative contain essentially the same information. Both 
together are expressions of the same disorder generated 
by the complex cooling procedure — which is subject to 
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the one chosen potential and its derivatives. It is the 
final positional disorder which matters, not the individ- 
ual contributions to the dynamical matrix. If the width 
of the lines between the particle positions shown in the 
snapshots (Fig. [I]) has been chosen proportional to the 
interaction force rather than the spring constant this was 
done mainly for artistic reasons. Resembling snapshots 
can be obtained, using as a scale for the line width either 
the spring constants or the trace rriij = M ix j x + M iy ^ y 
of the dynamical submatrices. Such pictures are direct 
visualizations of the dynamical matrix. 

The snapshots Fig. |l| reveal that the structure of the 
force network is very inhomogeneous, but isotropic on 
larger scale. They show also evidently that the quenched 
disorder is not much affected by the different symmetries 
(circle or square) and quench protocols. This corrobo- 
rates the robustness of the system properties with regard 
to the chosen protocol mentioned in the previous sec- 
tion. This is also supported by the different histograms 
regrouped in Fig. ||. 

Systems of the three protocols are included in first of 
both graphs Fig. ||(a) where we discuss the distribution of 
distances of interacting particles (inset) Tij/ro and of the 
tensile forces Uj (main figure) for relatively large systems 
containing N = 10, 000 particles. Comparison shows that 
there is essentially no dependence on quench protocol. 
Note that this is only a technical point and will make life 
easier when we analyse the eigenvalue spectrum of disk- 
shaped aggregates. The peak at r,j = tq in the inset 
corresponds to the peak for t t j = in the main figure of 
Fig. ||(a). Hence, most of the interactions have achieved 
to minimize their energy although this is not required 
for global mechanical stability. The distribution of the 
repulsive forces (negative tensions) is Gaussian. Inciden- 
tally, this is quite different from the recently repotted 
exponential force distributions in granular matterHH. 
Interestingly however, the tail of the force distribution 
becomes more and more exponential upon decreasing of 
the system size. We do not show this here, as the same 
effect is presented in the second panel (Fig. ||(b)) for the 
trace of the dynamical submatrices. Various system sizes 
as indicated in the figure are given here for the third pro- 
tocol at p = 0.925. This demonstrates that the tail of 
distribution which corresponds to very rigid contacts be- 
comes enhanced due to finite-size effect for systems below 
L w 20a. A very similar figure could also be given for 
the distribution of spring constants. Note that the peak 
on the left side of Fig. ||(b) is due to the slow variation of 
the spring constant for large distances and corresponds 
to distances seen in the shoulder on the right of the his- 
togram in the inset. 

In fig. H we show two of several attempts to analyse the 
correlations and fractal structure of force chains (or the 
spring constants) in view of extracting a characteristic 
size. We focus on strong repulsive forces and strong pos- 
itive spring constants. Only results from the first mea- 
surements are reported as both yield similar results. In 
a first step we obtain the network of all interactions with 



(negative) repulsive tensions above some threshold t u and 
compute then on these sets various correlation functions. 
The functions (cos 2 (t)) traced in Fig. ||(a) attempt to 
put the visual impression of linear force chains in quanti- 
tative terms. Here r denotes either the angle between the 
directors rii and rij of the links i and j (spheres) or the 
angle between the director rij of link i and the direction 
n-ij between the links i and j (squares). The difference 
between both methods is small. The correlation dies out 
only after about six oscillations. No significant difference 
have been found by increasing or decreasing the thresh- 
old t u . The decay of the envelope is exponential with a 
length scale of the order of the mean bead size. Hence, 
while rigid regions seems to be spatially correlated, this 
does not, apparently, introduce a new length scale- Note 
that the situation is similar in granular materialscj. 

The second part of Fig. ^ shows a direct attempt to elu- 
cidate the fractal structure of the network of quenched 
forces by means of the standard box counting technique. 
Here we count the number of square boxes of linear size 
b needed to cover all the links between beads carrying a 
repulsive force smaller than t u . For small b where every 
box contains only one link the differential fractal dimen- 
sion is zero. In the opposite limit where the boxes are 
much larger than the average distance between links, the 
differential fractal dimension must equal the spatial di- 
mension, i.e. df — 2. The power law slope df = 1 in 
between both limits indicates the typical distance where 
the contact network with tij < t u has a linear chain like 
structure. For example, we find a tangent with slope 
df = 1 at b ~ 6a for t u = —5. This distance is strongly 
increasing if we focus on more and more rigid subnet- 
works (decreasing t u ). Note, however, that there is no 
finite 6- window with df = 1 and in a strict sense there is 
again no characteristic length scale associated with lin- 
ear structures. A simple visual inspection of Fig. [l] that 
would suggest a characteristic length scale of the repul- 
sive force network much larger than the interatomic sep- 
aration is apparently incorrect and possibly caused by 
the natural tendency of our brains to emphasize linear 
patterns^. 

As a simple benchmark for the spatial correlation of 
the force network, we have distributed links randomly. 
The box counting of these uncorrelated links gives the 
dashed lines. Unfortunately, these are virtually identi- 
cal to the fractal dimension characterization of the force 
chains in the LJ systems and it is a very tiny difference 
(at least in this characterization) which is related to the 
spatial correlations. Obviously, this does not mean the 
forces are uniformly distributed as clearly shown by the 
snapshots and by the large number of oscillations in the 
correlation function. 

In summary, we provide evidence for a strong disper- 
sion in the dynamical matrix, i.e. in the local elastic 
properties of our systems. Although linear chain like 
structures exist apparently, this does not give rise to an 
additional characteristic length — defined as the screen- 
ing length in an exponential decay — much larger than 
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the mean particle distance. 



V. ELONGATION AND SHEAR: ELASTIC 
MODULI AND NON-AFFINE DISPLACEMENT 
FIELD 

A direct way of illustrating the failure of classical elas- 
ticity at small scales is to investigate the displacement 
field in a large, deformed sample. If the system is uni- 
formly strained at large scales, e.g. by compressing or 
shearing a rectangular simulation box, classical elasticity 
implies that the strain is uniform at all scales, so that 
the atomic displacement field should be affine with re- 
spect to the macroscopic box deformation. If this is true 
Eqns. ( ||) should provide reliable estimates of the Lame 
coefficients. Since the latter can be independently mea- 
sured in a computer experiment from the generated stress 
differencies using Hooke's law (Eqns. ( g)) this provides 
a direct and crucial test for the affinity assumption. 

Our numerical test proceeds in three steps: 

• Each initial configuration is given a first quick CG 
quench with quadruple precision to have a more 
precise estimate for the unstrained reference sys- 
tem. This is necessary since the stress differences 
we compute are very small for strains in the lin- 
ear response regime and the numerical accuracy of 
double precision simulation runs turns out to be 
insufficient. 

• Second, we perform an affine strain of order e of 
both the box shape and the particle coordinates. 
We consider both an elongation in x-direction 
(L x -> L x (l + e),r x -> r x (l + e)) and a pupaj 
shear using Lees-Edwards boundary conditionsES 
together with r x — > r x + r y e. (We refer here to 
the principal box particle coordinates (r x , r y ) of the 
periodic configurations.) 

• Finally, we quench with CG the configurations into 
the local minima while maintaining the strain at 
the boundaries. For given boundary conditions in 
the linear regime the solution must be unique. 

The differences of particle positions, forces and total en- 
ergies computed at each step are recorded. We stress that 
this procedure is technically not trivial and that great 
care is needed to measure physically sound properties. 

We have systematically varied e over several orders of 
magnitude from e = 1 down to e = 10~ 9 (not shown). 
The induced displacement field is reversible, linear in 
the amplitude of the imposed initial strain for an ini- 
tial strain in a window < e < 0.1 which decreases 
(for unknown reasons) somewhat with system size. Ob- 
viously, for higher strains the linearity of Eq. ( ||) must 
eventually break down. The departure from linearity be- 
low the given strain window is due to numerical accuracy. 



Using Eq. ( ||) for the stress tensor averaged on the 
whole system and the averaged strain field, we obtain the 
true macroscopic Lam coefficients A and /i. Within the 
given strain window the measured A and fi coefficients are 
strain independent, thus confirming the linearity of the 
elastic response in this range. The values of (x and A are 
presented in Fig. ^ where we have compared them with 
the affine field predictions. (Note that the Poisson ratio 
v w 2/3 is larger than 1/2 which is permissible in a 2D 
system as clarified in Sec. |l].) The coefficients relying on 
a negligible non-affinc field (open symbols, obtained from 
Eq. ( |9])) differ by a factor as large as two from the true 
ones. Clearly, a calculation taking into account the non- 
affine character of the displacements is necessary for dis- 
ordered systems. Indeed, it is simple to work out from the 
values given in the figures that, in a simple elongation, a 
finite energy fraction (A a + 2// )/(A+2/i) — 1 ~ 1/4 of the 
total strain can be recovered from the non-affine displace- 
ments. In a pure shear, the energy fraction fi a /fi ks 1/2 
is even larger; only the compressibility k = A + /i remains 
unchanged. Hence, in quantitative terms the non-affinity 
of the atomic displacements is not a negligible effect. 

The non-affine component Su of the atomic displace- 
ment field in large systems subject to an elongation in x 
direction is illustrated in the snapshot Fig. ||(a). A simi- 
lar snapshot holds in case of a pure shear (Fig. |(b)). In 
some regions, the displacement is much larger than ex- 
pected from a purely affine transformation (Note that 
even the non-affine part of the displacement field is 
strictly linear in e within the strain interval indicated 
above). Local displacements transverse to the direction 
of the elongation are allowed, and organize coherently 
into vortices. The transversal direction thus cannot be 
neglected, showing that the modelling approach put for- 
ward in Ref.li3 is not realistic. The crossover length £ 
mentioned in the Introduction manifests itself through 
correlated deviations from a purely affine displacement. 
Visual inspection tells us that the sizes of the vortices 
and £ w 30a are comparable. 

The two correlation functions presented in Fig. || con- 
firm this visual impression. The first one shows the cor- 
relation function C u (r) — (8u(r) ■ Su(0)). The striking 
anti-correlation for r <C 30a is in agreement with the size 
of the vortices seen in the displacement fields. That the 
displacement field is indeed correlated over a similar size 
is further elucidated in Fig. ||(b). Here we consider the 
systematic coarse-graining of the non-affine displacement 
field 

= jr. E f 11 ) 

of all Nj beads contained within the square volume el- 
ement Vj of linear size b. The mean-squared average 

U x (b) = (SUl// 2 is plotted versus the size of the coarse- 
graining volume element b. We have normalized the func- 
tion by its value at b = 1. The coarse-grained field 
decreases very weakly for b < 30a and only for much 
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larger volume elements we find the power law slope —1 
expected for uncorrelated events. As for symmetry rea- 
sons the total or mean non-affine field vanishes we have 
U x (b — > L) — > for very large volume elements. Apart 
from this trivial system size dependence U x approaches 
a system size independent envelope for b <C L, as can be 
clearly seen from the figure. 

Barely distinguishable functions have been obtained 
for U y (not depicted) which demonstrates the isotropy of 
the non-affine displacement fields which may also be in- 
ferred straight from the snapshots and appropriately cho- 
sen correlations functions. Similar characterizations can 
be obtained from a standard Fourier transform of 5u(x, y) 
and from the gradient fields defined on the coarse-grained 
field SU_y These are again not presented. 

In the rest of this section, we consider the local stresses 
generated by the applied macroscopic deformation. We 
show in Fig. the variance ((<J xy — cr yx ) 2 ) averaged on 
various boxes of size b. Here, the stress tensor o a p has 
been defined, as in classical mechanics, as the average 
force exerted in the a direction through the side perpen- 
dicular to the (3 direction of the volume element of size 
b. We see clearly in the Fig. that for a size b < 30a, 
this stress tensor is asymmetric, and that the asymme- 
try decreases exponentially to zero for larger sizes. In 
agreement with the discussion of section 0], this effect 
must be due to the lack of invariance of the energy un- 
der translations or rotations for volume elements of sizes 
< 30a. 

We note that the length scales observed in all these 
plots are relatively large compared to the particle size, 
but are not in the classical sense characteristic lengths ap- 
pearing in the exponential decay of a correlation function. 
An important consequence of the large spatial correla- 
tions is that calculations of Lame coefficients arc prone 
to finite-size effects for system sizes similar and below £. 
This explains qualitatively the peculiar system size de- 
pendence of the Lame coefficients reported in Fig. 0(b). 



VI. VIBRATION MODES 

In this section we discuss finally the eigenvalues and 
eigenvectors of the different systems we have generated. 
For each configuration the lowest (p < 1000) vibration 
eigenfrequencies and eigenvectors have been determined 
using the version of the Lanczos method implemented 
in the PARPACK numerical packageEJ. As stressed in 
the introduction we concentrate on the lowest end of the 
vibrational spectrum, since this is the part that corre- 
sponds to the largest wavelengths for the vibrations. 

We continue and finish first the discussion of disk- 
shaped aggregates which started with the snapshot 
Fig. ^(a). This being done we focus more extensively 
on the simpler periodic glassy systems (Protocol III) and 
discuss subsequently their eigenfrequencies and eigenvec- 
tors. This allows us to pay attention to central ques- 



tions of strongly disordered elastic materials without be- 
ing sidetracked by additional physics at cluster bound- 
aries (ellipticity, radial variation of material properties 
etc.). 



A. Eigenmodes of disk-shaped aggregates 

The first non-trivial eigenvalues (4 < p < 11) for 
the two protocols for disk-shaped clusters are shown in 
Fig. ^(a). The first three eigenmodes have vanishing 
eigenfrequencies because of two translational and one ro- 
tational invariance. Aggregates of different sizes are pre- 
sented as indicated in the figure. The frequencies are 
rescaled with the disk diameter 2R as suggested by di- 
mensional considerations or continuum theory. This scal- 
ing is roughly successful for all systems included. For 
the smaller systems (e.g., for the example with N = 732 
given) uj do not present the degeneracies of the continuum 
theory. If we increase the system size steps appear and 
the eigenfrequencies start to regroup in pairs of two fol- 
lowing the continuum prediction. This is well verified for 
the largest disk we have created containing N — 32, 768 
beads. The horizontal lines are comparisons with con- 
tinuum theory with appropriate density and where the 
Lame coefficients (and, hence, the sound velocities) have 
been taken accordingly from Table [II. The comparison 
of the two protocols (I and II) for N = 4096 shows that, 
perhaps surprisingly, the quench protocol does not mat- 
ter much if only the clusters are large enough and the 
excentricity sufficiently weak. Note that this is definitely 
not true for small disks where the excentricity matters. 

A schematic representation of the displacement fields 
for the lowest eigenmodes of a large disk-shaped clus- 
ter (R = 60) is given in Fig. [HJ For such large cluster 
the displacement fields essentially conform to standard 
elastic behavior. The two quantum numbers n and k 
are indicated for each mode. The snapshots demonstrate 
nicely the two-fold degeneracy for all the modes which 
are not axisymmetric (the first ten fields). We confirm 
that corresponding pairs are turned with respect to each 
other by an angle A9 = n/2k, e.g. for the first non- 
trivial modes p = 4 and p — 5 where k = 2 by an angle 
A9 = 45°. The disk wobbles in these modes between 
two perpendicular slightly elliptic shapes. The numerics 
returns arbitrary linear superpositions of eigenvectors as- 
sociated with the same eigenvalue. As we have at most a 
two-fold degeneracy this does, however, only amount to 
an arbitrary and common rotation of the pairs of snap- 
shots given. (The situation is more complex for periodic 
systems, see below.) Examples for the not degenerated 
axially symmetric case (k — 0) are the 'watch spring' 
p = 14 and the 'breathing mode' p — 15. The last mode 
is typically excited in Raman spectroscopyO. 

The next step is now to characterize the continuum 
approach as a function of system size. Our analysis 
presented in Fig. ^|(a) is similar to the finite-size stud- 
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ies of phase transitions and critical phenomena. The 
rescaled eigenfrequencies for given p are plotted versus 
the rescaled inverse system size £/2(R) in such a way that 
the vertical axis il 2 = {(ojR/ttc) 2 ) should become inde- 
pendent of the cluster properties (size, density, Lame co- 
efficients) in the limit of large systems. £ is the cross-over 
length defined in the previous section. We have taken 
£ = 30a. The horizontal lines correspond to the con- 
tinuum predictions for quantum numbers (n, k) as indi- 
cated. Both axes are dimensionless. For the p given in the 
figure only transversal modes are expected and, hence, we 
have used the transversal sound velocity c = ct (p) every- 
where. For p < 14 all modes should be two-fold degener- 
ated and one expects even (full symbols) and odd (open 
symbols) modes to regroup in pairs. This is born out 
for the larger systems; for 2(R) 3> £ the lowest frequen- 
cies even match quantitatively the predictions. There is 
no adjustable parameter left for the vertical scale. The 
crossover at £/2(i?) » 1 for the smallest eigenfrequen- 
cies justifies the (somewhat arbitrary) choice of the nu- 
merical value of £. Interestingly, the continuum limit is 
approached in a non-monotonous fashion and very small 
systems vibrate at higher frequencies. One cause for this 
is certainly the higher excentricity of smaller clusters. As 
we will see below, however, additional and more funda- 
mental physics plays also a role. 

Every data point corresponds to an average over an en- 
semble with identical operational parameters. The num- 
ber of configurations in every ensemble have been chosen 
such that the (not indicated) error bar is of the order of 
the symbol size. The dispersion of an individual mea- 
surement is much larger, however, for smaller systems 
and of the order of the frequency difference between sub- 
sequent modes ui(p + l) 2 — uj(p) 2 - As the problem seems 
to be strongly self-averaging, as one expects, the dis- 
persion between different representations of an ensem- 
ble goes strongly down with system size. Note that the 
diameters 2R and the densities p of each configuration 
in an ensemble vary somewhat and we have used in the 
averaging procedure the sound velocities associated with 
every specific sample density. This was done by means 
of interpolating the numerical values of the Lame coeffi- 
cients shown in Fig. |](a). The dispersion of sound veloc- 
ities within an ensemble is, however, relatively weak even 
though the Lame coefficients depend strongly on density. 

We have presented in the figure data from the first 
protocol of elaboration. Data from the second protocol 
looks quite similar. For small systems there arc differ- 
ences probably due to the higher eccentricity as men- 
tioned above. As some of the results presented here, such 
as the non-monotonous variation, could be due to spu- 
rious surface and eccentricity effects in circular systems, 
we now turn our attention for the rest of this section to 
the periodic glassy systems. 



B. Eigenfrequencies of periodic bulk systems 

Raw data for eigenfrequencies for systems generated 
following the third protocol are given in Fig. |](b) for two 
examples at p = 0.925. As there is no rotational invari- 
ance in a periodic box only the first two modes p — 1 and 
p = 2 vanish. The vibration frequencies do not display 
the degeneracies of the continuum in the smaller system 
with L = 32.9 (spheres). It appears that the finite-size ef- 
fects are much more pronounced in the eigenfrequencies 
compared to the weak effects discussed in Fig. |3|(b) on 
the stiffness. In contrast to small systems, the degener- 
acy steps are clearly visible for the largest configurations 
(square symbols) we have sampled with N = 40, 000. 
The quantitative agreement with continuum prediction 
is then satisfactory and deteriorates only slightly with 
increasing p, i.e. with decreasing wavelength X(p). Two 
configurations have been obtained in the latter case (open 
and full symbols). Interestingly, the self- averaging is such 
that both are barely different, even where they depart 
from the classical theory. 

Fig. ^|(b) shows the eigenfrequencies for bulk systems 
as a function of box size L in analogy to the characteri- 
zation presented above for disks. The horizontal axis is 
now£/L, the vertical axis £1 2 = ((ojL /2itct) 2 ) ■ While the 
continuum approach is somewhat smoother in the bulk 
case than for the disks essentially both sets of data shown 
in Fig. [9] carry the same message: They indicate that for 
system sizes L below £ the predictions of continuum elas- 
ticity become erroneous even for the smallest eigenmodes. 
The classical degeneracy of the vibration eigenfrequen- 
cies is lifted and the resulting density of states, becomes 
a continuous function. The approach of the elastic limit 
is again non-monotonous. This is not related to the dis- 
persion due to the discreteness of an atomic model, which 
would result into a monotonous approach to the elastic 
limit, as can be easily checked on one dimensional mod- 
els. As we have no surface effects here the effect must 
be due the frozen-in disorder. Obviously, the physics at 
play should also be relevant for the disk-shaped clusters. 

As a sideline we report here briefly that we have also in- 
vestigated the role of the quenched stresses on the eigen- 
frequency spectrum. This can be readily done by switch- 
ing off the contribution from the tensions in the dynam- 
ical matrix Eq. ( |^) . The corresponding artificial system 
appears to remain mechanically stable (i.e. all eigenval- 
ues remain positive), however, it does not correspond to 
any realistic interaction potential. A finite-size plot ana- 
log to the ones shown in fig. || has been computed for 
periodic systems. This yields a result qualitatively very 
similar to the curves presented here, albeit the crossover 
to continuum occurs for slightly smaller box sizes. The 
point we want to make here is two-fold: On the one side 
quenched forces matter if it comes to quantitative com- 
parisons and evaluation of a analytical model, on the 
other hand, they do not generate new physics. The role 
of quenched stresses is simply to maintain a local equilib- 
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rium in systems with strong positional disorder. The re- 
sults presented here, and particularly the deviation from 
classical continuum theory appears to be due to local 
disorder and simple harmonic but disordered couplings 
would give analogous results. 

As in the study of critical phenomena where finite- 
size effects reveal correlations at wavelengths X(L,p) in 
much larger systems where L>A one naturally expects 
that the results of Fig. ^ are also relevant for the de- 
scription of higher eigenmodes and their departure from 
continuum theory. We demonstrate in Fig. that the 
crossover to continuum is indeed characterized by the 
ratio of wavelength and the correlation length £. The 
numerical eigenfrcqucncies are rescaled by their ex- 
pectation n con t(n, m) from continuum theory (Eq. ( |^)) 
and plotted versus the inverse wavelength £/A cont (n, to) 
where the wavelength is inferred from Eq. ( |]) for the 
quantum numbers n and m associated with the mode in- 
dex p. Again we set (to some extend arbitrarily) £ = 30a. 
As can be seen, all the data sets obtained for various 
sizes collapse and confirm within numerical accuracy the 
choice of the scaling variables. Interestingly, two separate 
scaling functions appear for transverse and longitudinal 
modes and, for clarity, we have plotted both in two dif- 
ferent graphs. The crossover occurs at about X cont ~ £ 
for the transverse modes in agreement with the observa- 
tions in Fig. |9|(b). In contrast, about twice as large wave- 
lengths are required for longitudinal modes to obtain a 
satisfactory match with continuum theory as shown in 
Fig. 0(b). Both scaling curves are similar nevertheless 
and could indeed be brought to collapse by choosing a 
larger £ for the longitudinal waves. That £ depends some- 
what on the type of mode is not surprising. However, it 
would be of course more appealing if one would have a 
simple physical argument explaining the slower crossover 
for longitudinal waves. 



C. Eigenvectors: Noise and Correlations. 

Snapshots similar to those presented in fig. [To] could 
be presented for periodic bulk systems, however, as they 
turn out to be intricate linear superpositions of plane 
waves and as the existence of the continuum limit for the 
largest of our systems is now sufficiently demonstrated we 
focus in the reminder of this paragraph on the departure 
from the continuum prediction. 

In order to characterize the departure of the numer- 
ical eigenvector displacement fields v(p) from contin- 
uum theory we project them onto plane waves Ec ont (q), 
i.e. we compute their Fourier amplitudes A p (q) = 
(v(p)\v cont (q)). This is shown in Fig. for the eigen- 
vectors p = 3, 11 and 27. The average (...) is taken 
over the ensemble. The amplitudes are plotted versus 
q — p, i.e. we have shifted the abscissa axis horizontally 
in such a way as to emphasize the contribution of the 
p-th elastic mode to the computed eigenvector with the 



same number. Indeed, the main contribution is seen to 
be due to the plane wave (propagative) mode with the 
same mode number. The three particular eigenvectors 
considered in Fig. 0, belong to sets of four-fold degen- 
erate eigenstates. Hence, if noise could be discarded the 
projections onto plane waves would be of width four, cor- 
responding to an average over all possible random phases. 
Accordingly, the projection of eigenvectors belonging to 
an eight-fold degenerated set would have a width eight 
(not shown). As anticipated by our discussion of the 
eigenvalues the overlap between numerical and theoreti- 
cal eigenmodes deteriorates with increasing mode index 
p. This is seen from the decreasing peak height and the 
increasing width of the function A p (q— p) with increasing 
p. The enlargement of the peak to neighboring fre<™e&- 
cies suggests a scattering process in agreement witbliSO. 
Note the asymmetric character of the projections ampli- 
tudes, which must vanish for q < 3. Interesting, even for 
small p, the amplitudes do not completely vanish for large 
q — p, but become more or less constant. This indicates 
a Fourier transformed localized noise term, in agreement 
with the quasi-localized modes described in RefJiJ. As 
we shall elaborate now below, this is due to vortices in 
analogy to those depicted in Fig. ||. 

The next step consists in the construction of the noise 
field by substracting the contributions of the dominant 
peak from the numerical eigenvectors v(p) : 

Sv(p) = v(p) - ^ A p (q)v cont (q) (12) 
qev p 

where T> p is the set of 4 (or 8) plane waves q that con- 
tribute most to the Fourier decomposition of the mode p. 

Obviously, = \jYlq£v p A>(?) 2 should increase 

with the wavelength. We have computed the dimension- 
less ratio ||(?i>(jp )||/ ||i>(p)H and plotted this quantity in 
the inset of Fig. versus the wavelength X CO nt (p) ■ This 
curve is in qualitative-agreement with a scattering pro- 
cess of Rayleigh typeEH, where ||<5v(p)||/||v(p)|| oc A -2 . 
We have compared the data with the exponential decay 
cxp(— \ cont /30a). Interestingly, the characteristic wave- 
length defined here is equal to £. Our data do not allow 
to discriminate between both fits. The conclusion is thus 
that, whatever the origin of the noise (scattering process 
or not), the noise is small compared to the propagative 
theoretical mode when A 3> £. 

Let us now study the structure of the noisy part of 
the wavevector. Assuming a scattering process, it would 
be particularly interesting to determine the dependency 
of a possible mean free path into the wavelength X con t of 
the eigenmode. Two examples for eigenvector noise fields 
Sv(p) are presented in the figures g(c) and (d) for the 
modes p = 3 and p = 7 respectively. They are compared 
with the non-affine fields obtained for the same config- 
uration in an elongational and pure shear displacement 
field. The vortices are again the most striking features. 
The four fields given look indeed remarkably similar: The 
sizes and positions of the vortices are obviously highly 
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correlated. To put this in quantitative terms we consider 
correlation functions for the eigenvector noise fields de- 
signed in analogy to those discussed for the non-affinc 
fields. 

The correlation function C v (r) = (5v(r)\Sv(0)) is pre- 
sented in Fig. [[3] for two modes p = 3 and p = 11. As 
expected from the noise field snapshots we find again 
the anti-correlations similar to those presented for the 
non-affine fields in Fig. |](a). The anti-correlation ex- 
tends even to somewhat larger distances. We realize that 
the mode dependence while visible is weak. In the in- 
set we have plotted the first node of the anti-correlation 
(i(p). Also included is the wavelength corresponding to 
the mode number p. We find that (1 does not vary much 
with A, unlike the A 4 -dependence of the mean free path 
in scattering processes. Thus the noise displays a char- 
acteristic length (1 comparable to £ and independent on 
the mode p. We stress that the resulting participation 
ratiQ_of the noise is weak, in agreement with localiza- 
tiorit-i. However, the slow (non exponential) decay O^thc 
correlation function is in favor of delocalization as int-rEil 



VII. DISCUSSIONS. 

In summary, we have presented extensive simulations 
of mechanical and low-frequency vibrational properties 
of quenched amorphous disk-shaped aggregates and pe- 
riodic bulk systems. Two dimensional ensembles con- 
taining up to 40, 000 polydisperse Lennard- Jones parti- 
cles have been generated and analysed in terms of sam- 
ple size, density, sample symmetry and quench protocol. 
We have focused on systems with densities close to the 
zero-pressure state to have similar conditions for the two 
boundary symmetries studied. The eigenmodes of the 
structures are calculated by diagonalization of the dy- 
namical matrix and the eigenfrequencies are compared 
with the predictions from classical continuum theory 
where we concentrate on the low frequency end. The sec- 
ond key calculation we performed consists in macroscopic 
deformations (pure elongation and pure shear) of a peri- 
odic box in order to obtain the elastic moduli (Lame coef- 
ficients) and the microscopic displacement fields. These 
are in turn compared with the noisy part of the corre- 
sponding eigenvector fields. 

The central results of this study are as follows: 

1. The application of continuum elasticity theory is 
subject to strong limitations in amorphous solids, for 
system sizes below a length scale £ of typically 30 in- 
teratomic distances. This length scale is revealed in a 
systematic finite-size study of the eigenfrequencies and 
in the crossover scaling of modes p for fixed system size 
with the wavelength X(p). 

2. The success of the scaling of the eigenfrequencies 
with the wavelength demonstrates that £ is within nu- 
merical accuracy independent of the mode. The crossover 
behaviors of transverse and longitudinal modes are sim- 



ilar, although the latter is somewhat slower, i.e. slightly 
larger systems are required to match the same mode with 
continuum theory. 

3. The macroscopic deformation experiments demon- 
strate that the non-affine displacements of the atoms on 
local scale matter: a finite amount of energy is stored in 
the non-affine field and there is a large difference between 
the true Lame coefficients and those obtained from the 
dynamical matrix, assuming an affine displacement field 
on all scales. 

4. Below a length scale similar to £ both the non-affinc 
displacement field and the noisy part of the eigenvector 
fields displays vortex like structures. These structures are 
responsible for the striking anti-correlations in the vec- 
torial pair correlation functions of both types of fields. 
We have identified in this paper the non-affine field as 
the central microscopic feature that makes the contin- 
uum approach inappropriate. 

Inhomogeneities in local elastic coefficients, are an es- 
sential ingredient in sesietal recent calculations on dis- 
ordered elastic systemsEJo. Indeed, the origin of the 
departure from elastic behavior seen in the two key com- 
puter experiments is ultimately related to the local disor- 
der. This disorder is revealed by the structure of the force 
network frozen into the solid, as shown in Fig. Interest- 
ingly, weak finite-size effects arc evidenced in properties 
of the frozen local disorder like the pressure, the parti- 
cle energy and one of the Lame coefficients and in the 
histograms of forces, coupling constants and dynamical 
matrix trace. 

Surprisingly at first sight, we have been unable to iden- 
tify a sufficiently large length scale solely from the frozen 
local disorder. There must be a mechanism which ampli- 
fies these in such a way that that the deformation fields 
as the ones displayed in figure become non-affine on 
scales comparable to £ = 30a. Incidentally, this mecha- 
nism is not directly related to the mean free path that can 
be computed in diffusion processes^ as can be shown by 
the absence of any A-dependence. One marked difference 
is that .the mean free path is obtained in one-dimensional 
modelsE3'E£l, however the characteristic length is related 
in our study to vortices, implying a displacement in the 
transverse direction. A one-dimensional model of the di&. 
placement field in a two-dimensional medium as in Ref.Ea 
appears thus to be unrealistic. Unfortunately, we are not 
able at the moment to propose a definite relation between 
the size of the vortices and-the local properties of the sys- 
tem. Preliminary studiesEZl suggest a strong correlation 
between the vortices and the occurrence under mechan- 
ical sollicitation of nodes of stresses acting as bolts and 
forcing a displacement in the transverse direction. This 
must be related to the local-anisotropy of forces as al- 
ready been suggested in Ref.E3 in an other context. 

Interestingly, sizes-, similar to £, or somewhat smaller, 
are often invokedal!3, as typical of the heterogeneities 
that give rise to anomalies in the vibrational properties of 
disordered solids (glasses) in the Terahertz frequency do- 
main, the so-called boson peak. In particular, RefB con- 
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siders the existence of rigid domains separated by softer 
interfacial zones, not unlike those revealed by the non- 
afhne displacement pattern of Fig. ||(a). Our work offers 
a new vantage point on this feature. At the wavelength 
corresponding to these Thz vibrations, comparing the vi- 
brational density of states to that of a continuum, elastic 
model (the Debye model) is not necessarily meaningful. 

The present study documents the importance of sys- 
tematic finite-size characterizations for the computa- 
tional investigations of glassy systems well below the 
glass transition. It suggests that numerical investigation 
of vibrational properties should systematically make use 
of samples much larger than £, in order to avoid finite- 
size effects and to be statistically significant. Too small 
systems tend to be slightly more rigid and have higher 
pressures and system energies. 

Finally, let us mention that such vortices have been 
studied in disordered, mat erials in the context of large 
deformations and flowEJo. We show here that the same 
mechanism happens even in the elastic regime, for very 
small deformations in a quasi-static motion. Vortex like 
deformation patterns have also been identified, and asso- 
ciated with, in simulations of granular materialsEJ. Our 
studies clearly shows that very simple models, involving 
only conservative forces, are sufficient to reproduce such 
patterns. Disorder appears to be much more relevant 
than the 'granular' aspect (e.g. frictional terms) for this 
type of properties. 

The above conclusions are obviously subject to the 
conditions under which our simulations have been per- 
formed. The most serious limitation of this work is cer- 
tainly that we have only reported results on two dimen- 
sional samples. Of course, one expects correlations to be 
reduced in higher dimensions and finite-size effects should 
be less troublesome there. However, this study initially 
originated from an attempt to compute the vibrational 
modes of three dimensional clusters. Surprisingly, we 
have been unable to reach there the elastic limit even 
for systems containing 10, 000 atoms and had to switch 
to the simpler two dimensional case which is in terms of 
particle numbers less demanding even though the length 
scale £ might ultimately turn out to be smaller in three 
dimensions. Indeed, we believe that a systematic finite- 
size analysis of mechanical and vibrational properties in 
three dimensional amorphous bodies is highly warranted 
and we are currently pursuing simulations in this direc- 
tion. 

Other system parameters like the polydispersity index 
and specifically the quench protocol should be system- 
atically altered to further the understanding of the ori- 
gins of the demonstrated correlations. We have not seen 
within the accuracy of our data any systematic effect 
of the quench protocol on the the characteristic length 
£. However, more attention should certainly be paid to 
the ageing processes occurring in quenched samples as a 
function of system size. 
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-2.81 


24.8 


24.9 


31.9 


10000 
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-2.81 
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31.9 
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85 


76.9 
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0.16 
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24.1 


32.0 
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120 


109.9 


0.864 


0.11 


-2.80 


23.5 


23.1 


32.0 



TABLE I. Some properties characterizing disk-shaped clusters generated by a slow aggregation following protocol I with 
MD steps at T = 0.1,0.05,0.01,0.005 and 0.001. We have indicated the particle size N, the number of configurations in the 
ensemble M, the radius Ri of the initial sphere, the mean radius R of the final globule, the mean density p, the excentricity e 
obtained from the inertia tensor of the cluster, the interaction energy per particle E, the Lame coefficients A a and fi a obtained 
using eqns. ^ and the mean spring constant (rfj-Cjj). The mean tension (rijtij) and, hence, the pressure are too tiny to be 
measured accurately. 
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TABLE II. Some properties characterizing disk-shaped clusters generated with the fast quench protocol II which takes 
advantage of already quenched periodic bulk systems (protocol III). As in the first table we indicate the radius Ri of the initial 
disk, the number of configurations in the ensemble M, the (mean) particle number, the mean density p, the excentricity e, the 
interaction energy per particle E, the Lame coefficients X a and fi a and the mean spring constant (rfjCij). The final cluster 
radius is not given as it is essentially identical to Ri. Note that in the second protocol the particle number TV fluctuates (very 
weakly) around its mean value (TV) while it is a constant operational parameter in the first protocol. The excentricity is much 
smaller here than in the first more realistic protocol. The information contained in the last four columns is very similar to the 
one in the correponding columns of table I and table fnl despite the fact that different quench protocols have been used. 
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TABLE III. Some properties characterizing periodic bulk systems generated following protocol III with MD steps at 
T — 1.0, 0.5, 0.1, 0.05, 0.01, 0.005 and 0.001. In addition to properties also recorded in the previous tables we have included here 
the Lame factors A and obtained from a macroscopic deformation and the mean tension (rijtij). Note that {rfjCij) S> {rijUj). 
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FIG. 1. Representation of the network of quenched stresses in two small quenched Lennard- Jones particle systems in two 
dimensions: (a) a disk-shaped aggregate of diameter 2R w 32a containing TV = 732 particles (Protocol I) on the left and (b) 
a periodic bulk system with L = 32.9a and N — 1000 (Protocol III) on the right hand side. The line scale is proportional 
to the tension transmitted along the links between beads. The black lines indicate repulsive forces (negative tensions), while 
the red links represent tensile forces between the verticies. Both shown networks are very similar despite different symmetries 
and quench protocols. They are strongly inhomogeneous and resemble the pattern seen in granular materials. Zones of weak 
attractive links appear to be embedded within the strong skeleton of repulsive forces. 
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FIG. 2. Correlations of frozen in forces for large periodic systems: (a) pair correlation (cos 2 (t)) versus distance r. Here r 
denotes either the angle between the directors rii and rij of the links i and j (symbols) or the angle between the director m of 
link i and the direction mj between the links i and j (lines). The decay of the envelope is exponential with a length scale of 
the order of the mean bead size, hence, this does not introduce a new length scale. The dots correspond to an artificial force 
network obtained by shuffling randomly the tensions between existing links. (b) The number Nt of boxes of size b needed to 
cover the interactions transmitting repulsive tensions smaller than t u . The curves do not differ much from box counting results 
on sets of randomly drawn links (dashed lines). Different slopes are included for comparison. The slope -1 corresponds to linear 
chain like structures, the -2 slope to a compact structure in 2D. 
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FIG. 3. Histograms of contributions to the dynamical matrix: (a) comparison of the tension tij for all three protocols at 
N = 10,000. Inset: Histogram of distances rij/ro of interacting particles, (b) trace of matrix 

systems of various system sizes as indicated in the figure (protocol I). This demonstrates finite-size effects in the tail of the 
distributions of small systems with L < 20 corresponding to an enhancement of the skeleton of very rigid contacts. 
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FIG. 4. Lame coefficients A and fi obtained for periodic bulk systems: (a) density variation for TV = 10, 000, (b) size 
variation for p — 0.925. Full symbols correspond to the direct measurement using Hooke's law, open symbols are obtained 
using eqns. ^| Naturally, all Lame coefficients rise with density as the number and strength of pair interactions increases. Note 
that while /j, remains more or less constant A increases with decreasing L in a similar way as pressure and particle energy 
(tab. p3). 
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FIG. 5. Comparison of non-amne displacement field Su(r) with eigenvector field Sv (r) for periodic box of size L = 104a 
containing TV = 10, 000 particles: (a) nop-iaffine displacement field under elongation in x direction, (b) same for plain shear 
using Lees-Edwards boundary conditions^], (c) eigenvector field for p — 3 and (d) eigenvector field for p — 7. This confirms 
that the different noise fields are non-gaussian and are highly correlated in space and with respect to each other. Detailed 
inspection shows eddies like in turbulent flow. 
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FIG. 6. Characterization of non-affine displacement field <5u-field obtained by simple elongation: (a) correlations function 
C u (r) = {8u(r) ■ Su(0)}, (b) mean coarse-grained field U x (b)/U x {b = 1) versus the size b of the coarse-graining. Both functions 
become system size independent for large L. The first figure on the left side shows clearly an anticorrelation in agreement with 
the eddies seen in the snapshot fig. &a). 
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FIG. 7. Asymmetry of the stress tensor of the forces generated in simple shear (Lees-Edwards boundary conditions) for a 
system containing 10 000 particles. We measure the "microscopic stress" acting on a volume element as shown in the sketch 
on the right. In the main figure we plot the average mean-squared stress difference {{o xy — Oyx) 2 ) 1 ^ 2 versus the linear size b 
of the volume element. For small box sizes b < 30a we evidence power law behavior which crosses to an exponential decay at 
large volume sizes. 
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FIG. 8. Comparison of the first non-trivial eigenvalues for (a) disk-shaped clusters and (b) bulk systems confirming the 
predicted degeneracies for sufficiently large samples. In small systems, the degeneracy is lifted. The frequencies are rescaled 
with the system size (as indicated) and compared with the theoretical predictions (horizontal lines). The open symbols in (a) 
correspond to the slow quench (Protocol I), the full symbols to a rapid quench (Protocol II) showing that the results of both 
protocols become similar for large disks. The open and full squares in (b) correspond to two different configurations with the 
same parameters. 
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FIG. 9. Finite-size scaling of reduced eigenfrequencies Q.(p,L) — L)L/2tyct) 2 } versus £/L with £ = 30a : (a) 

disk-shaped aggregates from protocol I with L — 2R, (b) periodic bulk systems at p = 0.925 (protocol III). The mode index p 
increases from bottom to top, some of the p are specifically given (full symbols). In both cases the degeneracy is systematically 
lifted for small systems and the continuum prediction (given by the horizontal lines) is approached non-monotoneously. The 
pairs of quantum numbers (n, k) and (n, m) associated with the predictions are indicated in figure (a) and (b) respectively. 
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FIG. 10. First non-trivial eigenvectors (from p = 4 up to p — 15) of an aggregate of diameter 2R = 120a containing 
N = 10, 000 particles (protocol I). We have indicated the running index p and the pair of qantum numbers (n, k). The two-fold 
degeneracy for k > is obvious. Note that each pair is rotated by a polar angle A9 — -k /2k. The last mode represented is an 
axialsymmetric 'breathing mode' (k = 0). 
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FIG. 11. Scaling of (a) transverse and (b) longitudinal modes for different box sizes as indicated. The rescaled fre- 
quency 0,/Qcont = (uj 2 ) /uJcont is plotted versus the inverse wave length £/A cont . The theoretically expected wave length 
L/{n 2 + m 2 ) is given by the quantum numbers (n, m) associated with the mode index p. We have again set £ = 30a. 
The crossover to continuum theory occurs at A w £ for transverse modes and at about twice as large wave lengths for longitu- 
dinal modes. The success of both scaling plot suggests that £ is frequency independent for sufficiently large system sizes where 
A(p) <C L, but does depend the wave type. We regard this as the new central results of this paper. 



26 



N=10000 




FIG. 12. Construction of noisy eigenvector field for periodic box configurations with N — 10, 000, p = 0.925. Main figure: 
projection amplitude of empirical eigenvectors p = 3, 11 and 27 on the theoretical plane waves which are indexed with q with 
increasing frequency. Only the transverse modes are included for clarity. Insert: relative amplitude of noise as a function of 
wavelength A(p). The dotted line is a fit with exp(— \ CO nt/30a) in agreement with the estimation £ w 30a for the characteristic 
wave length. The long-dashed line is a fit with A -2 in agreement with a scattering process. 
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FIG. 13. Correlation function (5v (r) \5v (0)) of the eigenvector noise fields for p = 3 and 11 versus distance r. The 
curves are similar to the ones shown in fig. ^(b) and feature again prominent anti-correlations. Inset: Clip) characterizing the 
anti-correlation for N=10 000 (dotted line) and N=5 000 (bold line). Also included is the wavelength from Eq. ^ associated 
with each continuum mode (top line). 
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